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ABSTRACT 

A  theoretical  model  is  developed  for  the  wind  velocity 
field  in  the  atmospheric  near-water  layer  above  ocean  waves. 
A  turbulent  exchange  coefficient  due  to  the  wave-caused 
velocity  fluctuations  is  introduced  into  the  equations  to 
account  for  the  corresponding  turbulent  nature  of  the  motion 
and  the  resulting  energy  fluxes.   The  equations  of  motion 
for  air  are  derived  in  a  manner  to  allow  for  variation  of 
the  turbulent  exchange  coefficient  with  height.   A  stream 
function  is  introduced  and  a  numerical  solution  obtained  by 
utilization  of  the  Richtmyer   method.   The  profiles  of  wave- 
caused  Reynolds  stresses  and  their  resulting  energy  trans- 
fers are  discussed  and  compared  with  observational  data. 
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I.   INTRODUCTION 

Energy  and  momentum  transfer  to  water  waves  are  pri- 
marily dependent  upon  the  conditions  of  the  atmosphere 
immediately  adjacent  to  the  ocean  surface.   In  the  case 
of  air  flow  above  an  open  body  of  water  the  atmosphere- 
ocean  interactions  give  rise  to  wave  perturbations  in  the 
upper  layer  of  the  water  and  in  the  adjacent  atmosphere. 
It  has  been  suggested  that  the  wave  perturbation  charac- 
teristics be  used  to  determine  the  energy  and  momentum 
transfers  across  the  air-sea  boundary. 

Investigation  of  the  nature  of  these  energy  and  momen- 
tum transfers  have  been  the  purpose  of  numerous  studies 
in  recent  years.   Such  studies  of  wind-wave  coupling  were 
initiated  with  the  quasi-laminar  model  of  Miles  (1957), 
who  allowed  the  wave-induced  motion  to  interact  with  the 
shear  flow  only  at  the  critical  level.   Models  progressed 
through  the  studies  of  Phillips  (1957)  and  Benjamin  (1958) 
to  those  more  recently  of  Yefimov  (1970)  and  Reynolds 
(1972),  wherein  there  is  also  interaction  with  the  turbu- 
lent Reynolds  stresses. 

In  the  present  study,  wave  related  momentum  transfers 
were  investigated  through  the  use  of  a  numerical  model. 
This  study  is  similar  to  that  of  Yefimov  (1970) ,  but  it  is 
more  general  and  allows  for  dependence  of  the  mean  wind 
on  time  and  variation  of  the  turbulent  coefficient  with 
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height.   This  model  utilizes  an  initial  value  numerical 
procedure  similar  to  that  of  Newman  (1969). 

In  addition,  specific  data  derived  from  the  statisti- 
cal study  of  Davidson  and  Frank  (1973)  will  be  utilized 
to  define  solution  parameters  to  generate  the  numerically 
derived  stress  patterns. 
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II.   FORMULATION  OF  THE  MODEL 

The  basic  of  the  problem  is  the  equation  of  motion. 
Since  the  stress  is  expected  to  be  independent  of  y,  two- 
dimensional  motion  is  assumed  with  the  resulting  equation 
in  the  form: 

|^+  W  •  V\V=  -  i  VP  +  IF,  (2.1) 

where 


,2. 


Z     dx .  d  Z 

Here  x  is  the  horizontal  axis;  z  is  the  vertical  axis;  p 
is  the  density  of  air;  P  is  the  non-hydrostatic  pressure; 
u'  and  w'  are  the  turbulent  velocity  fluctuations.   The 
stationary  part  of  each  variable  (time  average)  will  be 
denoted  with  a  straight  line  above  it  and  the  nonstation- 
ary  part  (space  average)  by  a  waved  line.   The  stationary 
velocities  do  not  include  the  random  components  which 
would  disappear  with  a  suitable  time  average.   The  veloc- 
ity (W)  will  be  defined  as: 

W  =  (u  +  U)  a.  +   w  fk  (2.2) 

J 

3U 
where  U  =  U(z)  and  the  factor  -5—  will  be  neglected  in  the 

d  t 

development  of  the  equations.   For  simplicity  the  problem 
will  be  restricted  to  those  cases  of  a  single  wave  number 
This  neglects  the  effect  of  wave  self -interaction  which 
would  affect  other  waves. 
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The  turbulent  velocity  fluctuations,  u'  and  w'  are 
assumed  to  have  a  dual  nature.   They  are  generated  by  the 
gradient  of  the  mean  velocity  in  the  near  water  layer  and 
also  by  the  wave-caused  perturbations  of  the  velocity. 


The  tangential  Reynolds  stress  (u'w1)  will  be  expressed 

as : 

8U      ,du         9w, 


T  =  -u'w'  =  K  —■  +  V  (~  +   — )  ,       (2.3) 
dz  d  Z     dx 

where  u  and  w  are  the  horizontal  and  vertical  components 
of  wave-caused  velocity  perturbations.   The  factor  K  is 
the  coefficient  of  turbulent  exchange  which  is  related  to 
the  generation  of  turbulent  fluctuations  by  the  mean  veloc- 
ity gradient.   The  element  V  is  the  coefficient  of  turbu- 
lent exchange  which  is  related  to  the  wave-caused  velocity 

fluctuations.   It  will  be  assumed  that  V  >>  V   everywhere 

m 

except  in  the  thin  laminar  sublayer  in  which  we  shall  ne- 
glect the  effect  of  turbulent  viscosity.   Both  the  coeffi- 
cients, K  and  V,  will  be  assumed  to  be  constant  in  the 
x-direction  and  vis  a  function  of  height.   In  addition, 
the  coefficient  K  will  be  assumed  to  be  given  by: 

K  =  K  UA  z 

which  increases  linearly  with  height  and  leads  to  the  un- 
disturbed mean  velocity  profile 

U* 
U  -  —  In  (f-  +  1)  •  (2.4) 

o 
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Incorporating  the  continuity  equation  with  the  Bous- 

sinesq  approximation,  (V*W)  =  0,  the  equations  may  be 

written  in  component  form  as: 

3u  ,..,  „  ,  3U  1  3P  ,  3  ,  — t2n,  3  n/du*  ,„/9u  ,  3w.  ,  ,n  c. 
•7-r  +\V'  Vu+w-^—  =  -  —  -7T-  +  t— (-u'  )+  t-  K-r-  +v  (v-  +  -5—)  j  ,  (2.5) 
3t  dZ       p  3x     3x  3z    dZ      dz     dX 

dw  _m/  T7       1  3P  ,  3  ,„3U  ,   ,3u    3w.?,  3  ,  — r2v     ,  _  ,. 

■5—  +\V*  Vw  =  -  —  -r—  +  v-[K-r—  +  V  (-r—  +  -5—)  ]+  -r—  (-w'  ).  (2.6) 
3z  p  3z    3x  dz  3z    3x     dz 

The  vorticity  (£)  is  defined  to  be: 


3u    _    3\tf 
5    "    3z    "    3x    ' 


Normally    u 


-2    i    3   -r2 


w 


but  for  this  work  the  assumption  of: 


,2     ,2 
u    =  w 


will  be  incorporated  into  the  vorticity  equation  which  may 
now  be  written  in  the  form: 

15.  +   W.  U+   &    .Vu-  ^  -Vw  +  ^-(w-^)  =  ' 
3t        vt   3z       3x        3zv  3z; 


33u 


3z 


r„3U  ,  ,  3u  .  3w>,  ,  r 
-  [K5-  +v(t-  +  ^~)]-V[ 

2     dz      dz     dX        ^2 


3  3w 


3zz3z    3x 


(2.7) 


If  the  continuity  equation, 


3x    3z    U  ' 


is  used,  the  vorticity  equation  becomes 

2 


v—  +  U  -r—  +  w  +  u  -r—  +  w  -r^2- 

dt       3x       _  2       3x       3z 
3z* 


3U 


3u    3w. 


[KT7  +  ^TT  +  TT)]-V( 


33u 


33w 


'3z 


3z    3x' 


+  ^-)   (2.8) 


3zz  "~  "      3x23z    3x3 

The  products  of  the  perturbation  quantities  are  neglected 
except  where  they  affect  the  mean  flow.   This  is  consistent 
with  the  consideration  of  a  single  wave  number. 
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The  mean  quantities  are  balanced  separately  and  sub- 
tracted from  (2.8).   This  results  in  the  two-equation 
system : 


Ot         dX         r^_2  -N    2         °Z         OX  rs    2  r,  ~    3 


dz- 


3z- 


„  ir  +  w  ii  =  i 

ax      w  8z 


dZ 


dx^Bz    8x' 


2       dZ 


(2.9) 


(2.10) 


With  the  introduction  of  stream  function  notation,  the 
velocities  can  be  defined  as: 


u  = 


11 


w 


=  _  11 


9z  ■»  3x   ' 

which  transforms  (2.9)  and  (2.10)  into: 

I—  (l|)    +  i|j   )  =  -U(i|J     +  i|)    )  +  ip      U 
dt     zz     rxx  rXZZ      XXX      Tx 


(2.11) 


zz 


3z 


-    [v(lj>         -l(i       )]    -    v[ip  —    Tp  3     »     (2.12) 

2  rzz  rxx  irxxzz       rxxxx       ' 


ib   ib        +   ib    ib  -    ib    ib  -    ib   ib 

Z     XXX  Z     XZZ  X     ZXX  x     zzz 


Kibr 


(2.i3: 


z  zzz 

Yefimov  (1970)  simplified  the  solution  by  assuming  K  and  \ 
to  be  of  the  same  order.   This  assumption  will  be  incor- 
porated into  equation  (2.13)  by  replacing  K  by  V. 

The  analytic  solution  for  the  special  case,  U  =  0,  to 
(2,12)  and  (2.13)  has  been  considered  by  Yefimov  and 
Pososhkov  (1970).   Their  results  have  been  adapted  and  in- 
corporated as  a  test  of  the  numerical  technique  as  outlined 
in  Appendix  C. 
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The  stream  function  will  be  defined  as: 
r|»  =  tt»  (z,t)+  A(z,t)cos[k(x-ct) ]+  B(z,t)sin[k(x-ct) ]  (2.14) 

where  A(z,t)  and  B(z,t)  are  the  stream  function  amplitudes. 

Substitute  the  expressions  for  the  various  derivations 

of  ty    into  (2.12)  and  (2.13)  and  equate  coefficients  of  the 
cosine  terms  to  produce: 


f-(A    -  Ak2)  =  (U-c)tBk3  -  B   k]  +  BkU    + 
d  t    ZZ  zz  zz 


3z 


-    [v(A    +  Ak2)]  +  v[A   k2  +  Ak4]  . 
2       z  z  z  z 


(2.15) 


The  coefficients  of  the  sine  terms  are: 

f-(B    -  Bk2)  =  -(U-c)  [Ak3  -  A   k]  -  AkU    + 
3 t   zz  zz  zz 

—  [V(B    +  Bk2)]  +  V[Bzzk2  +  Bk4] 

The  mean  velocity  equation  is: 
k 


(2.16) 


IZ  z  zz 


(2.17) 


_[A  B    -A    B-AB+  AB    ]  =  V^ 
2l  z  zz     zzz      zz  z      zzz       o 

Equations  (2.15)  and  (2.16)  determine  the  wave-caused 
velocity  perturbations  and  therefore  are  the  basic  equa- 
tions for  solving  the  problem.   Equation  (2.17)  is  used 
for  the  calculation  of  the  mean  velocity  perturbation  u. 
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III.   SOLUTION  OF  THE  PROBLEM 

A.  INITIALIZATION 

To  initiate  the  solution  of  equations  (2.15)  and  (2.16), 
the  undisturbed  mean  velocity  was  assumed  to  be  in  the 

form  : 

U* 
U(z)  =  —  In  (|-  +  1)  ,  (3.1) 

o 

as  shown  in  Figure  1. 

The  functions  A(z,t)  and  B(z,t)  were  initialized 
through : 

A(z,0)  =  A(0,0)  sinh  [k (Hl-z) ] / sinh (kHl) ,         (3.2) 
and  B(z,0)  -  0  ,  (3.3) 

where  H  is  the  height  of  the  boundary  layer  and  HI  the 
height  one  grid  point  less  than  H. 

B.  BOUNDARY  CONDITIONS 

Yefimov  (1970)  defined  the  boundary  conditions  for  the 
lower  boundary  of  the  turbulent  layer  as: 

u(0,t)  =  akccos [k(x-ct) ] , 
w(0,t)  =  akcsin[k(x-ct) ] . 
In  terms  of  the  stream  function  the  problem  was  solved 
with  the  boundary  conditions: 

A(0,t)  =  ac,       A  (0,t)  =  -akc,  (3. A) 

B(0,t)  =  Bz(0,t)  =  0  ,  (3.5) 

A(H,t)  =  A  (H,t)  =  B(H,t)  =  B  (H,t)  =  0  .      (3.6) 
z  z 


17 


Z(rn) 


Figure  1.   Undisturbed  mean  velocity  profile  (m/sec). 
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Consequently,  the  boundary  layer  was  limited  to  a  certain 
height  H,  where  the  wave-caused  velocity  perturbations  were 

set  equal  to  zero. 

Yefimov  (1970)  studied  a  similar  problem  in  which  the 
upper  limit  of  the  boundary  layer  was  chosen  sufficiently 
high  so  that  it  was  reasonable  to  assume  the  conditions  of 
(3.6)  to  hold.   This  upper  limit  was  found  to  be  greater 
than  one  and  one-half  wavelengths  in  order  to  reach  a  sta- 
tionary solution. 

C.   RESULTANT  FACTORS 

As  a  consequence  of  the  initialization  functions  (3.1), 
(3.2)   and  (3.3);  the  boundary  conditions  (3.4),  (3.5)  and 
(3.6),  the  functions  A(z,t)  and  B(z,t)  were  time  stepped 
forward  until  a  steady  state  solution  was  reached.   At  each 
time  step  the  Reynolds  wave  stress  was  computed  from: 

T  -  -uw  f  -  |k  (BzA  -  AzB).  (3.7) 

According  to  equation  (2.17)  a  non-zero  Reynolds  stress  will 
result  in  a  mean  horizontal  velocity  (u)  of  the  wave-caused 
perturbations.   This  supplement  to  the  initial  undisturbed 

mean  velocity  may  be  computed  as: 

H 


u  =  - 

v 

o 


-  S  xdz  (3-8) 
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and  applied  as  a  correction  to  U  at  each  time  step  by: 

U*      z 
U(t)  =  —  In  (J-  +  1)  +  u  .  (3.9) 

o 


D.  FINITE  DIFFERENCING 

The  finite  differencing  was  accomplished  through  the 
use  of  centered  finite  difference  schemes  for  all  terms 
except  those  involving  friction.   These  terms  were  computed 
at  time  (t-At).   Examples  of  the  schemes  are: 

Az  =  \   AZ  (AK+1  "  AK-1>  '  (3-10> 

and 

Azz=  77  <AK+1  "  2AK  +  Vl)  '         (3'11) 
Az^ 

As  a  result  of  the  finite  difference  schemes  and  the 
boundary  conditions,  the  variables  at  the  second  grid  point 
are  specified  by: 

B(Az  ,  t)  =0  . 
In  addition  the  difference  schemes  and  the  boundary  condi- 
tions specify  that  the  stream  function  amplitudes  A(z,t) 
and  B(z,t)  remain  constant  at  both  the  bottom  and  top  two 
grid  points . 

E.  NUMERICAL  SOLUTION 

Equations  (2.15)  and  (2.16)  were  solved  utilizing  the 
principle  of  Richtmyer  (1957)  with  a  forward  time  step  to 
start  the  marching  procedure.   To  avoid  solution  separation, 
the  problem  was  restarted  every  50  time  steps  with  a  finite 
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difference  scheme  developed  by  Matsuno  (1966).   This  was  a 
two-step  iteration  to  simulate  the  Euler  backward  differ- 
ence method. 

All  calculations  were  performed  to  double  precision 
accuracy  and  continued  until  the  Reynolds  stress  (3.7) 
reached  a  steady  state  condition  which  is  defined  by: 


/  H  2\l/2 

=   I   E   [T(z,t)  -  T(2,t-1)]/  1        < 


k  -  I  \  •    ii  (z,t  )    -  t(z,L-  l»  ]-  !      :    l.  x  10  5 


o 

A  space  step  of  0.25  meters  and  a  time  step  of  0.0125 
seconds  were  utilized.   Additional  accuracy  was  attained 
in  some  cases  by  reducing  the  space  step  to  0.125  meters 
with  a  corresponding  time  step  of  0.00312  seconds. 
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IV.   COMPUTATIONAL  PROCEDURES 

A.   VARIABLES 

The  forecast  formulas,  as  stated  in  (2.15)  and  (2.16), 
allow  for  five  independent  variables.   These  are  wave 
amplitude  (a),  wave  number (k)  ,  dynamic  velocity  (UA)  ,  rough- 
ness length  (z  ),  and  the  coefficient  of  turbulent  exchange 
o 

(V). 

Yefimov  (1970)  assumed  that  the  slope  of  the  waves  was 
small,  so  that  the  expressions  for  the  velocities  on  the 
boundary  predicted  by  the  theory  of  waves  of  small  slope 
could  be  applied.   This  assumption  of  a  small  slope  is 
applied  by  establishing  the  relationship  ak  =  .01,  which  is 
incorporated  in  the  model  as  a  means  of  computing  the 
initial  wave  amplitude  A(0,0). 

The  resulting  variables  k,  U^,  zq,  and  V  were  then  con- 
sidered systematically  by  first  determining  the  range  of 
variation  of  each  of  the  variables.  Then  central  values 
were  selected  and  used  as  example  cases  by  holding  three 
variables  fixed  and  computing  test  cases  of  the  remaining 
variable  over  its  full  range. 

Kraus  (1967)  studied  the  relationship  of  the  roughness 

length  to  that  of  the  mean  wind  where  roughness  length  is 

given  by: 

z   =  —  U. 
o    G   * 
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The  constant,  m,  was  determined  by  different  researchers 
with  values  ranging  from  0.0125  to  0.08.   A  central  value 
of  m  =  0.042  was  selected  and  the  roughness  length  was  com- 
puted from: 


z   =  0.00428  u\   , 
o  * 


(4.1) 


for  the  case  where  z   was  dependent  upon  U..   This  reduces 

the  independent  variables  to  l'A,  k  and  V. 

The  majority  of  the  experiments  were  conducted  with 

z   =  0.0001  meters  after  Hill  (1962). 
o 

B.   MODEL  CONSIDERATIONS 

Although  the  value  of  the  coefficient  of  turbulent  ex- 
change is  unknown,  it  is  reasonable  to  assume  that  it  can 
be  variable  with  height  as  discussed  by  Yef imov  (1970) . 
This  possibility  has  been  incorporated  in  the  model  and 
several  cases  have  been  investigated  where  V  was  allowed 
to  vary  with  height. 
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V.   RESULTS 


A.   EXPERIMENTAL  RESULTS 

The  number  of  variables  and  their  large  range  of  varia- 
tion necessitated  numerous  numerical  experiments.   To  re- 
duce the  number  of  experiments  as  much  as  possible,  the 
variables  were  restricted  to  a  limited  range.   The  selection 
of  a  variation  increment  was  guided  by  the  necessity  to 
choose  a  sufficiently  large  increment  while  remaining  with- 
in the  scope  of  full  investigation  of  the  problem.   The 
variables  and  the  range  of  values  used  are  shown  in  Table  1. 
These  were  used  in  the  various  combinations  as  shown  on 
Figures  2  -  10. 


TABLE  I.   Range  of  the  variables  considered  in  the 
numerical  experiments. 


kCm"1) 

0.050 
0.075 
0.100 
0.125 
0.150 
0.175 
*0.200 
0;225 
0.250 
0.275 
0.300 


U^ (m/sec) 

0.010 
*0.020 

0.040 

G.060 

0.080 
*0.100 

0.150 
*0.200 

0.250 

0.300 

0.350 

0.400 


V (m  /sec) 

1.100 
*0.240 
0.200 
0.110 
0.020 
0.011 
0.002 


z^  (m) 
— o 


0.10 
0.01 
0.001 
0.0005 
*0.0001 
0.00005 
0.00001 
0.000005 
0.000001 


*central  values 
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The  existence  of  a  critical  level  was  an  important 
factor  in  the  determination  of  the  profile  of  the  solution, 
The  critical  level  is  defined  as  that  height  where  phase 
speed  is  equal  to  the  mean  wind  speed.   Since  the  phase 
speed  is  inversely  proportional  to  the  wave  number(c  = 


/G/k) ,  the  phase  speed  will  decrease  as  wave  number  in- 
creases.  Therefore,  for  a  fixed  mean  wind  profile,  an 
increase  in  wave  number  will  decrease  the  height  of  the 
critical  level  (if  existing) .   In  the  cases  of  fixed  wave 
numbers,  an  increase  in  UA ,  or  decrease  of  zq  should  de- 
crease the  height  of  the  critical  level. 

The  first  case  considered  was  that  of  a  variable  wave 
number.   In  this  case  the  results  were  greatly  affected 
by  the  value  chosen  for  UA  and  therefore  two  experiments 
were  conducted.   In  the  first  experiment  (Figure  2)  the 
wave  number  was  varied  with  11*  =  0.20  m/sec.   The  second 
experiment  (Figure  3)  was  identical  but  UA  =  0,020  m/sec. 

Some  noteworthy  features  of  these  two  tests  are: 

a.  the  decrease  in  magnitude  and  variation  of  T  as 
U.  was  decreased; 

b.  the  decrease  in  magnitude  and  variation  of  X  as 
k  was  increased; 

c.  a  maximum  value  of  x  occurs  at  approximately  one 
meter  above  the  wave  surface; 

d.  the  increase  in  T  at  a  certain  wave  number  when 
UA  =  0.20  m/sec.   It  is  significant  that  this  was  the 
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first  case  in  which  the  critical  level  was  within  the  layer 

where  x  varies  significantly. 

The  second  case  considered  was  that  of  different  mean 

wind  profiles.   The  mean  wind  given  by  (3.1)  is  dependent 

upon  by  U.  and  z  .   Since  U.  and  z   are  considered  to  be 
r  J       *  o  *  o 

interrelated  there  are  numerous  ways  to  vary  the  mean  wind 
profile.   In  this  study  three  experiments  were  conducted. 
The  first  case  (Figure  4)  was  that  of  a  fixed  wave  number, 
roughness  length  and  turbulent  coefficient.   In  the  second 
case  (Figure  5) the  wave  number  and  turbulent  coefficient 
remained  fixed  but  the  roughness  length  was  allowed  to  vary 
according  to  (4.1).   Hence,  both  dynamic  velocity  and  rough- 
ness length  varied  in  this  case.   The  third  case  (Figure  6) 
was  that  of  fixed  dynamic  velocity,"  wave  number  and  turbu- 
lent coefficient,  while  roughness  length  was  allowed  to 
vary. 

Some  notable  features  are: 

a.  the  increase  in  magnitude  of  the  wave-related 

stress  as  z.  was  decreased; 
o 

b.  the  increase  in  magnitude  of  the  wave-related 
stress  as  U^  was  increased  until  U   is  approximately  equal 
to  0.25  m/sec.   This   is  the  point  where  the  critical  level 
is  within  the  layer  where  X    varies  significantly. 

The  third  case  to  be  considered  was  that  of  various 
turbulent  coefficients  (Figure  7)  while  the  other  variables 
were  held  f  ixed . 
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A  noteworthy  feature  is  the  increase  in  magnitude  of  T 

2 
as  V  decreased  to  0.020  m  /sec.   This  is  followed  by  a 

decrease  in  the  magnitude  of  T  as  V  continued  to  decrease. 

In  the  fourth  case  the  turbulent  coefficient  was 
allowed  to  vary  with  height  (Figure  8) .   In  the  first  ex- 
periment the  turbulent  coefficient  was  decreased  with 
height,  to  10  meters,  according  to  the  relationship: 
V(z)  =  0.5  sinh  (H-z) /sinh (H) . 

On  the  subsequent  experiments  V  increased  linearly  with 

2 
height  (Figure  8)  until  V  =  0.24  m  /sec  and  then  remained 

constant  thereafter.   The  results  of  the  first  experiment 
are  shown  in  Figure  9  and  those  of  the  increasing  V  experi- 
ments are  shown  in  Figure  10.   A  similar  case,  but  with  V 
held  constant,  is  shown  in  Figure  10  for  comparison. 
Some  notable  features  are: 

a.  the  similarity  in  the  profiles  of  the  cases  of 
increasing  V  with  that  of  decreasing  V; 

b.  the  presence  of  a  maximum  at  approximately  0.5 
meters  above  the  wave  surface; 

c.  the  rapid  change  from  positive  to  negative  stress; 

d.  the  decrease  in  magnitude  of  X    as  V  increases  to 
greater  heights. 

B,   COMPARISONS  WITH  DATA 

Additional  experiments  were  conducted  to  compare  the 
results  of  this  model  with  that  of  the  data  analyzed  by 


2? 
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0.003280 

0.179 

43.50 

23 

0.000079 

0.273 

29.00 

23 

0.000130 

0.273 

none 

Davidson  and  Frank  (1973).   The  initial  input  parameters 

(Table  II)  were  combined  with  a  constant  turbulent  coeffi- 

2 
cient  (v  =  0.240  m  /sec)  and  the  calculations  time  stepped 

forward  until  a  steady  state  was  reached.   The  resulting 

stress  profiles  and  amplitudes  are  shown  in  Figures  11-13. 


TABLE  II.   Dynamic  velocity,  wavelength  and  roughness 
length  values  determined  by  Davidson  and 
Frank  (1973)  with  corresponding  wave  numbers 
and  critical  levels. 

-     critical 
Case    U. (m/sec)   A (m)    z  (m)    k(m   )    level(m) 

1.  0.200 

2.  0.167 

3.  0.121 


Appendix  E  outlines  the  procedure  for  obtaining  data  from 
the  reported  results  of  Davidson  and  Frank  (1973)  and  also 
the  procedures  for  interpretation  of  the  numerical  results. 
The  elements  for  comparison  have  been  displayed  in 
Table  III. 

A  significant  feature  of  Table  III  is  the  high-degree 
of  correlation  of  the  magnitude  of  the  velocities  obtained 
from  the  variance  spectra,  phase  amplitude  and  numerical 
amplitude  procedures.   An  equally  significant  feature  is 
the  correlation  of  the  magnitude  of  the  stresses  obtained 
by  the  velocity  cospectra  with  those  from  the  numerical 
calculations.   It  should  be  noted  that  the  sign  of  the  co- 
spectra  stress  is  dependent  upon  the  wave  phase.   The 
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differences  in  the  signs  should  therefore  not  be  considered 
as  a  significant  discrepancy. 

The  comparisons  are  the  first  of  this  kind  for  data 
from  the  natural  regime.   The  agreement  is  significant. 
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Figure  2 


Distribution  of  the  wave-caused  stress  (x  in 
m2/sec2)  with  height  for  U^  =  0.2  m/sec,  V  = 
0.240  m2sec2,  z   =  0.0001m  and  k  equal,  respec- 
tively, to:     °    1)  0.050,  2)  0.075,  3)  0.100, 
A)  0.125,  5)  0.150,  6)  0.175,  7)  0.200, 
8)  0.225,  9)  0.250,  10)  0.275,  11)  0.300.   The 
arrows  indicate  the  height  of  the  critical  level 
and  the  identifying  wave  number. 
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Figure  3.   Distribution  of  the  wave-caused  stress  (t  in 
m2sec2)  with  height  for  UA  =  0.02  m/sec,  V  = 
0.240  m2/sec,  z   =  O.OOOlra  and  k  equal,  respec- 
tively, to:     °    1)  0.050,  2)  0.075,  3)  0.100, 
A)  0.125,  5)  0.150,  6)  0.175,  7)  0.200, 
8)  0.225,  9)  0.250,  10)  0.275,  11)  0.300. 
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Figure  4.   Distribution  of  the  wave-caused  stress  (t  in 
m2/sec2)  with  height  for  k  =  0.200m~l,  V  = 
0.240  m2/sec,  z  =  0.0001m,  and  U*  equal,  respec- 
tively, to:     °   1)  0,01,  2)  0.02,  3)  0.04, 
4)  0.06,  5)  0.08,  6)  0.10,  7)  0.15,  8)  0.20, 
9)  0.25,  ]0)  0.30,  11)  0.35,  12)  0.40.   The 
arrows  indicate  the  height  of  the  critical  level 
and  the  identifying  dynamic  velocity. 
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Figure  5.   Distribution  of  the  wave-caused  stress  (t  in 

m2/sec2)  with  height  for  k  =  0.20m-1,  V  =  0.240 
m^/sec,  z  =  0.00428  U^2  and  U^  equal,  respective- 
ly, to:   °   1)  0.01,  2)  0.02,  3)  0.04,  4)  0.06, 
5)  0.08,  6)  0.10,  7)  0.15,  8)  0.20,  9)  0.25, 
10)  0.30,  11)  0.35,  12)  0.40.   The  arrows 
indicate  the  height  of  the  critical  level  and 
the  identifying  dynamic  velocity. 
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Figure  6.   Distribution  of  the  wave-caused  stress  (x  in 
m2/sec2)  with  height  for  U   =  0.100  ra/sec, 
k  =  0.200  m~l,  V  =  0.2A0m2/sec  and  z   equal, 
respectively,  to:   1)  0.1,  2)  0.01,   3)  0.001, 
4)  0.0005,  5)  0.0001,  6)  0.00005,  7)  0.00001, 
8)  0.000005,  9)  0.000001. 
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Figure  7.   Distribution  of  the  wave-caused  stress  (t  in 
m2/sec)  with  height  for  k  =  0.200  n"1,  U*  = 
0.100  m/sec,  z   =  0.0001m,  and  V  equal, 
respectively,  °to:   1)  1.100,  2)  0.200, 
3)  0.110,  A)  0.020,  5)  0.011,  6)  0.002. 
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Figure  8.   Variation  of  the  turbulent  coefficient  with 
height,  respectively,  to:  1)  decreasing, 
2)  increasing  to  5m,  3)  increasing  to  10m, 
4)  increasing  to  15m,  5)  increasing  to  20m, 
6)  increasing  to  25m. 
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Figure  9.   Distribution  of  the  wave-caused  stress  (x  in 
m^/sec^)  with  height  for  U .  =  0.100  m/sec, 
k  =  0.200  m"1,  z   =  0.0001m  and  V  decreasing 
with  height . 
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Figure  10.   Distribution  of  the  wave-caused  stress  (t  in 
m^/sec^)  with  height  for  U  =  O.lOOm/sec,  k  = 
0.200m--'-,  z  =  0.0001m  and  V  increasing  with 
height,  respectively,  to:  1)  5m,  2)  10m,  3) 
15m,  4)  20m,  5)  25m.-  Curve  6)  is  the  case  of 


constant  V  with  the  same  U 


k  and  z   values 
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_  Figure  11 . 


Distribution  of  the  wave-caused  stress  (t  in 
m2/sec2)  with  height  for  k  =  0.179m~l,  UA  = 
0.200  m/sec,  z  =  0.003280m,  and  V  =  0.240m2/ 

2  ° 

sec  . 
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Figure  12.   Distribution  of  the  wave-caused  stress  (t  in 
m2/sec2)  with  height  for  k  =  0.273m  *- ,    U^  - 
0.167m/sec,  z  =  0.000079m,  and  V 
.2  ° 


0.240m2/ 


sec1 
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Figure  13.   Distribution  of  the  wave-caused  stress  (t  in 
m2/sec2)  with  height  for  k  =  0.273m-1,  U. 


0.121m/sec,     z 
.2  ° 


=    0.000130m,     and    v=0.240    m* / 


sec 
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VI.   CONCLUSIONS 

The  Reynolds  wave  stress  of  the  perturbation-velocity 
field  in  the  atmospheric  near  water  layer  above  waves  has 
been  obtained  from  a  numerical  model.   The  resultant  Rey- 
nolds stress  can  be  used  to  compute  the  energy  which  is 
transferred  by  the  average  motion  of  the  air  to  the  wave 
fluctuations  and  vice  versa.   Thus  a  positive  stress  indi- 
cates the  wind  energy  maintains  the  wave  fluctuations  sig- 
nifying a  transfer  of  energy  from  the  wind  to  the  waves. 

The  maximum  Reynolds  stress  and  greatest  energy  trans- 
fer was  found  to  exist  for  the  largest  mean  wind  profile 

(largest  U.  and  smallest  z  ),  the  lowest  wave  numbers  and  a 
*  o 

2 
turbulent  coefficient  of  approximately  0.020  m  sec.   In 

addition,  the  maximum  stress  occurred  at  approximately  one 

meter  above  the  wave  surface. 

A  striking  result  is  the  increase  in  the  magnitude  of 
the  Reynolds  stress  as  the  mean  wind  velocity  approaches 
that  of  the  wave  phase  speed  in  the  lowest  levels  of  the 
boundary  layer.   For  those  cases  where  the  critical  level 
occurs  in  the  region  of  the  maximum  stress  fluctuations 
there  is  a  significant  increase  in  the  magnitude  of  the 
stress.   Thus  the  energy  is  greatest  when  wave  velocity  is 
close  to  the  mean  wind  velocity  in  the  lower  levels. 

Another  significant  result  is  the  change  which  occurred 
in  the  stress  profile  when  the  turbulent  coefficient  was 
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allowed  to  vary  with  height.   The  net  effect  of  this  change 
is  a  decrease  in  the  total  energy  transfer. 

The  comparison  of  numerically  generated  velocities  and 
stresses  with  those  of  measured  data  yielded  a  very  high 
degree  of  correlation.   These  results  will  justify  future 
numerical  research  and  data  collection  to  provide  new  in- 
sights into  the  air-sea  energy  transfers.   Specifically,  the 
correlation  of  a  numerical  model  with  measured  wind  data 
could  determine  the  characteristics  of  the  turbulent  coeffi- 
cient due  to  the  wave-caused  motion.   Concurrent  studies 
being  conducted  by  Aurand  (1973)  and  Strieker  (1973)  may  be 
considered  as  sources  for  additional  data  to  correlate  the 
statistically  and  numerically  derived  stress  patterns. 
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APPENDIX  A 
BAROTROPIC  INSTABILITY  TEST 

Neuman  (1969)  considered  a  simple  barotropic  numerical 
model.   The  object  of  the  test  was  to  determine  the  growth 
rate  as  a  function  of  the  wave  number.  The  similarity  be- 
tween the  present  study  and  that  of  Neuman  allowed  for  a 
check  of  the  present  model  for  the  special  case  of  no 
f rict  ion . 

However,  to  accomplish  the  check,  several  modifications 

were  made  to  the  present  model.   These  are: 

A(z,0)  =  cos  („  )  ,  and 
n 

V(z)  =  0. 

The  mean  wind  velocity  profile  is  given  by: 

U(z)  =  -U   tanh  (y/y  )  (Al) 

as  shown  in  figure  14. 

The  procedure  required  marching  the  calculations  for- 
ward in  time  until  the  growth  rate  stabilized.   The  point 
of  maximum  Reynolds  stress  for  three  different  times  were 
then  compared.   The  magnitude  of  the  stress  at  these  points 
was  found  to  be  growing  exponentially,  as  was  expected. 

As  an  additional  check,  the  growth  rate  of  the  Reynolds 

stress  was  evaluated  and  compared  to  the  growth  rate  as 

given  by  the  function: 

2U  n 

exp( — )t  .  (A2) 

o 
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Figure    14.       Mean    wind    profile    for    U    - 
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In  this  case  the  stress  at  the  maximum  point  was  equated 

to  (A2)  and  the  equation  solved  for  n  using: 

U   =9.0  m/sec  , 
o 

y   =  4 . 0m  . 
'  o 

The  resultant  growth  rate  was  0.181.   This  is  more 
accurate  than  that  of  Neuman  who  found  a  growth  rate  of 
0.160,  but  it  is  less  than  that  of  Betchov  and  Criminale 
(1967)  who  set  the  rate  at  n  =  0.19.  The  latter  difference 
is  because  the  present  solution  has  been  solved  with  a 
finite  boundary. 
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APPENDIX  B 
EXAMINATION  OF  THE  VISCOSITY  TERMS 

In  this  test  the  initial  field  of  the  cosine  stream 
function  amplitudes  was  computed  utilizing  equation  (3.2). 
At  this  point  the  mean  wind  and  phase  speed  were  set  to 
zero.   Then  the  forecast  equations  were  marched  forward 
until  a  steady  state  was  reached. 

These  numerical  values  were  then  compared  with  the 
analytical  values  derived  from  the  solution  of  the  equa- 
tion : 


A      +  2A    +  A  =  0  .  (Bl) 

zzzz      zz 


The  boundary  conditions  utilized  in  the  solution  of  (Bl) 
are  similar  to  those  given  by  (3.4)  and  (3.6).   These 
conditions  are: 


A(0)  =  2c  ,  A  (0)  =  -2kc  , 


A(H)  =  A  (H)  =  0  . 
z 


The  solution  of  (Bl)  is  given  by 


A(z)  =  c,  sin(kz)  +  c„  cos(kz)  +  c_  kz  sin(kz) 

+  c,  kz  cos(kz)  ,  (B2) 
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where : 


ci  -  i  iV0)  -  ca  k]  • 


c2  =  A  , 


C3  ''      ~    Hk  sin  (Hk)  l"z 


[A  (0)sin(Hk)+  A(0)cos(Hk)  + 


c,  sin(Hk)+  c,  Hk  cos  (Hk) ]  , 


)  + 


A„(0) 


Hk 


-]sin(Hk) 


c   .  *  5|B  W    ([A(0 
sin  (Hk)-H  k  \ 

+  [^^-  +  A(0)k  cot(Hk)]  cos(Hk)])  . 

The  comparison  of  the  numerical  results  with  the 
analytical  results  derived  from  (B2)  indicated  that  the 
minimum  occurred  at  the  same  height  for  both  cases.  How- 
ever, the  magnitude  of  the  numerical  minimum  was  smaller 
than  that  of  the  analytical  results. 
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APPENDIX  C 
EXAMINATION  OF  THE  CASE  OF  NO  MEAN  WIND 

Yefimov  and  Pososhkov   (1970)  derived  an  analytic  solu- 
tion to  the  forecast  equations  (2.15)  and  (2.16)  for  the 
special  case  U  =  0.   This  solution  is  given  by: 


-2zk  .    -zk//2a      ,zk  .        /„,\ 
+  ae         cos  ( )  ,     (CI) 

/la 


where: 

Vk 
a  =  — 
c 

The  condition  U  =  0  was  established  and  the  numerical 
solutions  of  (2.15)  and  (2.16)  were  computed  to  a  steady 
state  condition.   The  results  of  the  analytical  and  numeri- 
cal calculations  of  stress  were  compared  for  several  dif- 
ferent wave  numbers  and  for  different  values  of  the  turbu- 
lent exchange  coefficient. 

The  results  of  the  stress  obtained  by  both  methods 
were  very  similar.   In  all  cases  compared  the  analytical 
stress  reached  a  minimum  at  a  lower  height  and  somewhat 
smaller  magnitude  than  the  numerically  derived  stress. 
This  discrepancy,  however,  has  been  primarily  attributed 
to  the  scale  of  the  grid  increment. 
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APPENDIX  D 
EXAMINATION  OF  THE  RICHTMYER  SCHEME 

This  test  was  devised  solely  to  test  the  method  of 
solving  equations  as  given  by  Richtmyer  (1957).   The  test 
was  conducted  by  equating  the  right  hand  side  terms  of 

7T  Z 

equation  (2.15)  to  sin  ( )  where  w  =  (H-3)Az.   The  result 

w 

is : 

(7T  "  k  }  8T  '  sin  (w~)  • 

a  z 
and  upon  solution  yields: 


A  =    -rr  91 2~  At  •  (D1) 

[(£)   +  k  ] 
w 


The  numerical  solution  was  marched  forward  one  complete 
time  step.   Then  the  minimum  value  of  the  stream  function 
amplitude,  A(z,t),  was  compared  with  the  value  derived 
from  (Dl) .   The  result  was  a  comparison  of  the  values 
to  within  three  percent  of  the  difference. 
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APPENDIX  E 
EVALUATION  OF  ANALYZED  DATA 

Davidson  and  Frank  (1973)  reported  observed  values  of 
wave  heights  and  velocity  fluctuations  at  two  levels  above 
the  waves.   The  results  from  their  spectral  and  joint  pro- 
bability density  function,  conditional  mean  function  (JPDF- 
CMF)  analyses  are  shown  in  Figures  15  and  16. 

The  velocity  variance  spectra  graphs  show  plots  of 
n  •  cf>  (n)  vs.  ln(n)  where  n  is  the  frequency  of  the  wave  in 
Hz.   The  area  under  the  significant  spikes  of  the  curve 
(shaded  areas  of  Figure  15)  can  be  used  to  determine  the 

value  of  the  square  of  the  wave-caused  velocity  fluctua- 

2       2 
tion  (u*   or  w'  ).   These  areas  may  be  computed  by: 

Area  =  n  •  <j>(n)  *  A  ln(n)  ,  (El) 

where  A  ln(n)  is  defined  at  the  significant  peak. 

The  velocity  cospectra  graphs  (Figure  15)  are  similar 
to  those  above  but  the  area  in  this  case  is  indicative  of 


the  wave-caused  Reynolds  stress  (u'w').   These  results 
have  been  presented  in  Table  III  for  comparison. 

Davidson  and  Frank  (1973)  also  presented  phase-ampli- 
tude results  for  the  same  data  as  above.   These  results 
(Figure  16),  which  may  be  interpreted  directly,  are  shown 
in  Table  III  for  comparison. 
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The  numerical  maximum  amplitude  has  been  computed 
from  the  steady  state  values  of  the  stream  function  ampli- 
tude by: 

9         9   1/9 

(E2) 


u 


'  -  (A  2  +  B  2)1/2 

2         Z 


and 


w'  =  k(A2  +  B2)1/2 


(E3) 

which  follow  as  a  result  from  equations  (2.11)  and  (2.14). 
The  numerical  results  were  computed  with  the  small  wave 
slope  assumption  where  ak  =  0.1.   However,  this  exact  con- 
dition is  not  normally  met  in  nature.   Therefore,  the 
values  obtained  from  equations  (E2)  and  (E3) ,  as  well  as 
the  numerically  derived  Reynolds  stress,  must  be  scaled 
to  reflect  the  actual  wave  slope. 

The  calculation  of  a  scaling  factor  first  involved  the 
determination  of  the  actual  amplitude  (a  ) .   This  may  be 
interpreted  directly  from  the  "waves"  diagram  in  Figure  16 
The  amplitude  used  for  the  numerical  calculations  may  be 
computed  by: 


a„  = 


0.1 


N     k 

The  resulting  velocity  scaling  factor  (Fs)  may  now  be  ex- 
pressed by: 


Fs  = 


(E4) 


and  the  quadratic  stress  scaling  factor  (Fq)  by 

2 


Fq 


(E5) 


'N 
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